Nodal Protectorate: A Unified Theory of the afr-plane and c-axis Penetration Depths 

of Underdoped cuprates 
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We formulate a model describing the doping (x) and temperature (T) dependence of the afc-plane 
and c-axis penetration depth of a cuprate superconductor. The model incorporates the suppression 
of the superfluid density with underdoping as the system approaches the Mott-Hubbard insulating 
state by augmenting a d-wave BCS model with a phenomenological charge renormalization factor 
that is vanishingly small for states away from the nodes of the d-wave pair potential but close to 
unity in the vicinity of the nodes. The c-axis penetration depth is captured within a model of 
incoherent electron tunneling between the CuCh planes. Application of this model to the recent 
experimental data on the high-purity single crystals of YBa2Cu30e+i implies existence of a "nodal 
protectorate", a fc-space region in the vicinity of the nodes whose size decreases in proportion to x, 
in which d-wave quasiparticles remain sharp even as the system teeters on the brink of becoming an 
insulator. The superfluid density, which is extremely small for these samples, also appears to come 
exclusively from these protected nodal regions. 



I. INTRODUCTION 

It is believed that to understand the high-T c cuprate 
superconductors, one must understand the problem of 
doping a Mott insulator i Conversely, understanding the 
manner in which superconducting order gives way to 
the antiferromagnetic insulator at half filling would pro- 
vide important clues to the solution of this fundamen- 
tal challengei2i2i^ Indeed, those cuprates which have 
the lowest concentrations x of dopants (so-called "under- 
doped" cuprates) possess many of the most enigmatic and 
mysterious experimental properties. In particular, as x is 
reduced, the zero-temperature pair-potential maximum 
Ao grows while the zero-temperature superfluid stiffness 
p s (0) and transition temperature T c decreased In addi- 
tion, the underdoped cuprates possess a "pseudogap" 7 
in the single-particle excitation spectrum that persists 
above T c . Taken together, these suggest that the way 
superconductivity is destroyed as x — > is very unusual. 
Unfortunately, there has been a paucity of data on such 
very underdoped cuprates, in part due to sample prepa- 
ration difficulties. In addition, such materials are often 
very disordered, complicating analysis of their intrinsic 
physical properties. Recently, however, Liang et al. 8 have 
devised a way to vary x in the doping phase diagram for 
a single crystal of YB^CuaOe+s (YBCO) in an essen- 
tially continuous fashion. Briefly, their method utilizes 
the fact that the effective doping on the Cu02 planes 
is governed by both the oxygen concentration and the 
degree of oxygen ordering in the CuO chains. By vary- 
ing the latter variable (via room temperature annealing) , 
these authors can alter the effective doping (and thus T c ) 
to explore the doping dependence of physical quantities 
(such as the c-axis penetration depthp^tiS) in the strongly 
underdoped regime x — ► 0. 

Motivated by the recent c-axis penetration results of 
Hosseini et al., 9 - 10 in this Paper we develop a model 
aimed at capturing the overall phenomenology of the dop- 



ing and temperature dependence of the afr-plane and o 
axis penetration depths (X a b and A c , respectively). Our 
model combines incoherent tunneling in the c-direction 
with a phenomenological momentum-dependent charge 
renormalization factor inspired by the idea of Ioffe and 
Millie 1 1 to account for anomalous doping dependences. 
On a more fundamental level we are interested in under- 
standing the constraints that this data imposes on the 
models of strongly correlated electron matter describing 
the d-wave superconductor on the verge of becoming a 
Mott-Hubbard insulator. Our results indicate that the 
appropriate effective theory must exhibit a "nodal pro- 
tectorate" consisting of regions in the vicinity of the d- 
wave gap nodes where the quasiparticles remain well de- 
fined even as x — > 0. 

To begin, let us review the overall phenomenology of 
the penetration depth in the cuprates. For the moment, 
we are interested in general trends (especially in the 
strongly underdoped regime x —> 0) in the penetration 
depth. To the extent that such data is available the 
a6-plane penetration depth X a b(x, T) exhibits the follow- 
ing behavior: 
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where we have expressed \ a b(x,T) in terms of the as- 
sociated superfluid stiffness p^ b (x,T) measuring the free 
energy cost to a variation of the order-parameter phase 
in the ab plane. Here, d denotes the distance between 
copper-oxygen layers. For YBCO the relevant param- 
eters are a ~ 244meV and b ~ 3.0i 12 i 14 We shall as- 
sume that this phenomenology holds for all underdoped 
cupratesi£*i&. The coefficient of T in Eq. (|TJ is nearly 
doping independent, in the sense that although it varies 
as the gap amplitude Ao varies, it does not vary as 
strongly as p^ b (x,0) i.e., linear in x. The T-linear term 
is well known and is understood to be due to the ex- 
citation of quasiparticles out of the condensate near the 



2 



nodes of the d-wave order parameter. 17 The fact that this 
relation holds to very low dopings (along with other ex- 
periments sensitive to nodal physics, see, e.g., Ref. [Isj ) 
indicates that near the nodes elementary excitations in 
cuprates are well described by conventional BCS quasi- 
particles even for x — > 0. 

The doping dependence of p® h (x,T), on the other 
hand, represents one of the central mysteries in the 
cuprates^ The linear in x reduction of p° b (a;,0) as one 
approaches half filling must be attributed to the prox- 
imity of the Mott-Hubbard insulating phase. However, 
most theoretical treatments (such as the RVB-typeiS and 
Gutzwiller projection22i2i approaches) that can account 
for the observed x-dependence of p^ b (x,0) also predict 
a strongly x-dependent prefactor b in Eq. JJJ, in dis- 
agreement with the data. Thus, the central theoreti- 
cal problem appears to lie in constructing a model that 
would make only a small fraction ~ x of all electrons par- 
ticipate in the superconducting condensate while at the 
same time preserve the simple BCS character of the nodal 
quasiparticles. In the following we resolve this problem 
by postulating that, at least for the purposes of study- 
ing the superfluid stiffness, the underdoped cuprates can 
be described by the BCS theory augmented with a phe- 
nomenological constraint that only Cooper pairs com- 
posed of electrons with momenta in the vicinity of the 
nodal points effectively couple to the external electro- 
magnetic field. We implement this idea by extending the 
"effective charge renormalization" concept introduced in 
this context by Ioffe and Millis. 11 

This model, described in a more detail below, is de- 
signed to reproduce the afr-plane phenomenology embod- 
ied in Eq. QJ. What makes us believe that it might 
be of more general validity is the fact that it also re- 
produces the doping dependence of the c-axis superfluid 
stiffness, p c s (x,T), once we adopt a model for interlayer 
coupling that gives the correct (nonlinear) temperature 
dependence of this quantity. The phenomenology along 
the c-axis is tantalizingly similar to the ab plane case and 
can be written in the following way»2ii2i22i2£ 



ing to experimental data, the parameters determining 
this factor are extracted. In Sec. 11111 we consider a tun- 
neling model of c-axis transport and augment the result- 
ing expression for A c with the same charge renormaliza- 
tion factor in an analagous way. It is shown that the T 
and x dependences of our expression for A c agree quali- 
tatively (having made certain assumptions) with Eq. J2Jl . 
In Sec.JIIII we fit our expression for A c to the results of 
Refs.l MlCl) . In Sec. El we make some concluding remarks. 
Certain calculational details are relegated to Appendices 
E|amd El 



II. afe-PLANE PROPERTIES 

The starting point of our calculation of the a6-plane 
penetration depth A a &(T) is the following Hamiltonian 
for a d-wave superconductor coupled to an applied in- 
plane electromagnetic field: 



H — ifpair + -Hint , 
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p c s (x, T) cx \~ 2 (x, T) ~ Ax a - BT a , 



(2) 



where -ff pa ir is the standard BCS pairing Hamiltonian 
with the single-particle energy ej. = — 2t(cos k x a + 
cos k y a) — p with a the lattice spacing. The pair po- 
tential Ak = |Ao(cosfe x a — coskya); we shall always 
assume that the maximum pair potential Ao is approx- 
imately temperature independent. Here, Hi nt refers to 
additional interactions not captured by i? pa ir which de- 
scribe the physics of the proximity to the Mott-Hubbard 
insulating phase at half filling. Our strategy will be to 
compute A~ b 2 by first neglecting H m t and computing the 
usual BCS result for the superfluid density in a d-wave 
superconductor. Then, we shall include the residual in- 
teractions H mt via a phenomenological charge renormal- 
ization motivated by the work of Millis et alSL and Ioffe 
and Millisiii The calculation of the penetration depth 
for -ffp a i r is standard and is presented in Appendix 1X1 for 
completeness. The final result may be expressed as 



where the exponent a is, within the experimental accu- 
racy, the same for both x and T and close to 2, while 
A, B are again x and T independent constants. The 
fact that the temperature dependence is nearly quadratic 
(as opposed to linear) points to the incoherent coupling 
between the copper-oxygen planes, as discussed by pre- 
vious authorsi22i2£ In the following we shall clarify the 
specific conditions under which such a nearly quadratic 
T-dependence arises for the incoherent tunneling model. 
More importantly we show that within our scheme for 
the effective charge the same power law also necessarily 
governs the doping dependence of p c s {x, 0), in agreement 
with Eq. ©. 

This Paper is organized as follows: In Sec.[Hj we aug- 
ment the standard BCS expression for A a b with a phe- 
nomenological charge renormalization factor. By appeal- 
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tanh -f3E k , 
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with E'k = \/ 1\ + A k the BCS <i-wave excitation spec- 
trum and (3 = 1/T the inverse temperature. The factor 
Z£ will be discussed shortly. 

We shall take Eq. (j^J to be the starting point of our 
phenomenological analysis of the a6-plane penetration 
depth data. It contains a sum over all momentum vec- 
tors k in the Brillouin zone and has been cast into a 
form where the zero temperature value A o fc(0)~ 2 and the 
finite temperature correction SX a f,(T) = A^ fc 2 (0) — A~ fc 2 (T) 
are treated on an equal footing; in particular they are 
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FIG. 1: Schematic plot of the assumed form for Z^, showing the "nodal protectorate" region (shading) of the Brillouin zone 
where states contribute to the formation of the condensate for a cuprate at a particular doping x. The black lines are the 
constant energy contours in the Brillouin zone (which do not vary with doping). Near optimal doping [x = 20), electrons in a 
large region around the node contribute to the Meissner response. As the doping is reduced this region is progressively reduced, 
leaving a small "patch" near the nodes where superconductivity remains robust. We remark that the c-axis penetration depth 
measurements of Ref. lj] were performed on extremely underdoped samples with effective dopings x that are approximately 
represented by the leftmost panel. 



both dominated by the regions in the fc-space close to 
the nodes of the d-wetve gap function. As discussed in 
Appendix ^ for — 1 this expression is fully equiva- 
lent to well-known expressions for X 2 b that have appeared 
previously in the literature 

Following Ref. we have introduced in Eq. Q a 
phenomenological charge renormalization factor Z^ by 
means of the replacement 



(5) 



This is done to incorporate the effects of the interactions 
contained in H ln t responsible for the gradual demise of 
superconducting order as the Mott-Hubbard insulating 
phase is approached near the half filling. The main dif- 
ference between our approach and that of Ref. |n| is that 
we have incorporated Z^ into the full expression for \ 2 ab , 
as opposed to the temperature dependent part only. This 
is in keeping with our philosophy of treating both com- 
ponents on the equal footing. At this stage we do not 
attempt to justify the inclusion of from microscopic 
considerations; we merely state that such charge renor- 
malization is not prohibited by any general principle (and 
is known to occur e.g. in the Fermi liquid theory). We 
offer some discussion of the possible origins of Z^ in Sec. 
IV. 

In the absence of a microscopic theory for Z\^ we must 
rely on experimental data to infer the behavior of this 
quantity. We begin by reiterating that, in weak coupling 
BCS theory Z\^ — 1 and a direct evaluation of Eq. Q 
yields a result which does not conform to the observed 
afe-plane phenomenologjsii 

In particular BCS theory yields a correct linear-T de- 
pendence of SX a b(T) but the wrong doping dependence 
of A~ b 2 (0), which would scale with the total number of 
electrons (1 — x) in disagreement with Eq. To rec- 
tify this discrepancy the form of Z^ must be modified. 
Since at low temperatures the T-dependence of 5X a b(T) 
comes from thermally excited quasiparticles near the four 



nodal points of a d-wave gap, it is clear that in order to 
preserve this correct temperature dependence Z^ must 
remain close to unity in the nodal regions. To suppress 
the overall magnitude of A~ b 2 (0) it then follows that 
must be small outside the nodal regions. In addition, to 
conform to the A~ b 2 (0) ~ x scaling the size of the "patch" 
in which « 1 must scale with x. 

The a6-plane phenomenology thus dictates the follow- 
ing form for the charge renormalization factor 



Z forE k <E c 
for E k > E c 



(6) 



with Zq a constant of order 1 and E c ~ x the charac- 
teristic cutoff energy that is chosen to obtain the correct 
magnitude of the zero-temperature superfluid stiffness. 
As illustrated in Fig. 2J this choice of Z^ effectively re- 
stricts the sum over the entire Brillouin zone in Eq. Q 
to the immediate vicinity of the nodes. At low temper- 
atures, such a choice affects the diamagnetic response 
(which is sensitive to the entire Fermi surface) but not 
the paramagnetic response which governs the tempera- 
ture dependence. 

To see how this reproduces the assumed phenomenol- 
ogy, we proceed by inserting the assumed form for Z^ 
into Eq. We make an assumption, which we verify 
momentarily, that E c is sufficiently small (i.e. E c -C Ao) 
for underdoped cuprates that it is permissible to linearize 
ek and Ak in the vicinity of the four nodesi^ We imple- 
ment this linearization by introducing a local coordinate 
system (fci, £2) centered at the nodal point with axes ro- 
tated 45° with respect to (k x , k y ). Expanding to leading 
order, we have £k — > vpk\, Ak — > Ua&2 with Vf, «a be- 
ing the Fermi and gap velocities at the node, respectively. 
We thus obtain 
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where E^ — y/vpk\ + v\k\ is the linearized "Dirac" 
spectrum of quasiparticles. It is now useful to rescale 
the integration variables as vpk\ — > ki and t>A&2 — ► &2 
and pass to polar coordinates. The angular integral is 
trivial and we obtain (inserting h and fee) 
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2TtTi 2 d va 
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tanh sech — 

2 2 2 



{E c -ik B T]n2) 



(8b) 



where the second line applies in the regime T -C -Be- 
lt is clear that if we take E c oc x, then this repro- 
duces Eq. (JTJ. To make a rough estimate of how E c 
must vary with the T c of an underdoped sample for this 
picture to apply, we simply assume that T c may be de- 
termined by the T at which A~ h 2 (T) in Eq. l|8b")l equals 0. 
This criterion, which was also used by Lee and WenjM 
gives E c ~ 2.8fc B T c = 0.24T c meV/K. Comparison with 
the experimental data shown in Fig. [2] illustrates that 
this overestimates T c by about a factor of 2.3 due to the 
fact that the actual data deviates significantly from the 
straight line at higher temperatures. To account for this 
deviation, and for future reference, we revise our estimate 
to read 



E c ~ 6.4fc B T c = 0.55T c meV/K. 



(9) 



This implies that even at optimal doping, T c ~ 90K, 
the energy scale E c is of the order of maximum gap 
Ao ~ 45meV. This validates our assumption that for 
underdoped YBCO the linearized approximation should 
lead to quantitatively correct results. Similarly one can 
estimate the doping dependence of E c by comparing Eqs. 
J3J and JSEJl. For YBCO this yields 



E r 



2na — Z, 



o 



[219.0meV]a 



(10) 



where we have taken Zq(vf/va) = 7, a value relevant 
to this material pU* In Sec. IIVI we shall see that a similar 
value of E c also provides the best fit to c-axis penetration 
depth data. 

There is an important caveat regarding the preced- 
ing analysis due to the fact that it relies on on the ap- 
proximate formula Eq. (|8b(l . Indeed, the exact formula 
Eq. H8a|) does not possess a T c in the sense of Eq. Il8bll : 
One can easily see that A~ fc 2 (T) > for any T. This is be- 
cause we have taken A k to be nonzero and temperature- 
independent as implied for underdoped cuprates by var- 
ious experiments such as ARPES^i or tunneling^ For 
non-zero pair potential within mean-field theory, the con- 
densate is depleted only through thermal excitation of 
quasiparticles; this leaves a small fraction of electrons 
in the condensate at arbitrary temperatures. In other 
words within mean-field BCS theory the only way to kill 
the condensate is to drive A^ to zero. In Fig. ©, we plot 
Eq. I|8a|) (solid line) and Eq. (|8b|l (dashed line) as a func- 
tion of the normalized temperature T/E c . One can show 



— Theory Equation (8a) 
■ ■ Low T Extrapolation 
♦ Experimental Data 




FIG. 2: Plot of normalized superfluid stiffness 
(Xab(0)/\ab(T)) 2 with \ ab (T) given by Eq. (g£j) (solid 
line) and Eq. 18bt (dashed line), respectively, as a function 
of the normalized temperature T/E c . Although they agree 
well at low temperatures, the deviation of Eq. 18bl from 
Eq. I8al becomes significant at high temperatures. As 
noted in the text, Eq. 18at never reaches zero, indicating 
that the quasiparticle-excitation mechanism of depleting the 
condensate can never fully destroy superconductivity if the 
pair potential remains nonzero. Diamonds represent the 
afo-plane data (YBCO 6.6) of Ref. |T3l scaled in such a way 
that the low-T linear part agrees with the theoretical model. 



that the solid line approaches zero only asymptotically as 
1/T 3 . In any real system, once the superfluid stiffness be- 
comes sufficiently small due to quasiparticle excitations, 
fluctuation effects (e.g., phase fluctuations^) will rapidly 
destroy the condensate. Thus, it is reasonable to assume 
that the relationship © holds approximately for real sys- 
tems; i.e., it provides an estimate for how the physical T c 
depends on the zero-temperature superfluid stiffness (pa- 
rameterized by E c ). In the next section, we compute the 
c-axis superfluid stiffness within a tunneling model and 
will find similar behavior. 



Finally we may ask whether there is another form of 
Zk that might agree with the experimental data. One 
possibility is that, for some reason, Z^ does not van- 
ish outside the nodal patch but goes to a small value 
proportional to y/x. This would clearly reproduce the 
observed afr-plane doping dependence^ We shall see be- 
low, however, that such a form would not produce the 
correct doping dependence for A c ; the latter appears to 
dictate that Z^ vanishes outside the patch. Also, we have 
only considered situations when Z^ depends on the mo- 
mentum through the energy E^. While this assumption 
appears very natural one can envision scenarios in which 
Zk would depend on the individual components of k. 
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III. c-AXIS PROPERTIES 

In the present section, we apply the philosophy of 
Sec. ^ to the problem of the c-axis penetration depth, 
for which new data in the strongly underdoped regime 
has recently been obtainedAiS 

A. Tunneling Hamiltonian 

To model the electronic transport in the c-direction, 
we adopt the following tunneling Hamiltonian: 

where c r:m a annihilates an electron at position r in layer 
m with spin projection a. The matrix element t T con- 
nects states in adjacent layers and is distinct from the 
quantity t in Sec. [H] The total system Hamiltonian, 
then, consists of a sum of intralayer Hamiltonians given 
by Eq. © coupled by H T - 

To compute the c-axis superfiuid stiffness, we derive 
via the Kubo formula the response to an AC electric field 
E(r, m, t) = E(r, m)e tuJt in the c-direction. To do this, 
we recall that such an electric field may be incorporated 
via the Peierls substitution 

t r t r e i ^ A(r - m) , (12) 

where A(r, m) — —E(r, m)e iuJt and d is the interlayer 
spacing. The c-axis current density j(r, m) is given by 

j(r,m) = -iet r (cJ i7n+1)CT c rimi<T - h.c.) 

e 2 d t 
+— t r A(r, m){cl m+l a c Tim ^ + h.c). (13) 

The two terms in Eq. I|13|) have the familiar form of the 
paramagnetic and diamagnetic contributions to the cur- 
rent, respectively. We proceed to compute the associated 
contributions to the conductivity to leading order in per- 
turbation theory in the matrix element i r . We shall fur- 
thermore assume that the in-plane Green functions are 
given by the usual rf-wave BCS forms: 

( r rC_ Pimii (T) Cpm , iT (0)) = <y tn , ro /F(p,T), (14a) 
(TVc PiroiCT (r)c pimV (0)) = -5 m , m ,G(p,T), (14b) 

with G(p,u>) and F(p,u>) given in Eq. (|A12(1 and 
Eq. i|A13(l . respectively. 

Of principal importance is our choice for the form of 
the tunneling matrix element t r . Allowing for nontrivial 
r-dependence is essential, as a spatially uniform t T (so 
that in-plane momentum is conserved during the tun- 
neling process) yields a T-linear correction to the pen- 
etration depth at low temperatures. This may be un- 
derstood by noting that such purely coherent tunneling 
between layers simply yields a three-dimensional d-wave 



superconductor. As noted previously^S^ the absence of 
such linear behavior points towards incoherence in the 
c-axis tunneling. Before proceeding, we remark that an- 
other possibility to obtain non-linear-T dependence is 
the proposal of Xiang and Wheatleyi2& They find that 
by including a momentum dependence to the tunneling 
matrix element (arising from the structure of the cop- 
per and oxygen orbitals involved in tunneling), one finds 
a T 5 dependence in disagreement with the overall ob- 
served temperature dependence [Eq. It is possible 
that such a T-dependent contribution is present but sim- 
ply masked by the dominant T 2 behavior; however, we 
shall not include this possibility in the following. 

Our strategy for finding a form for t p _k that gives the 
experimentally observed T-dependence is to assume that 
interplane disorder scatters the momentum states tun- 
neling from layer to layer. We assume that the disorder- 
averaged matrix element £ k = 0. For the disorder- 
averaged product of tunneling matrix elements we choose 

t*J^ = (2TT) 2 S(q)T k 2 , (15a) 

r k 2 = ^ k2/K \ (15b) 

with t j_ being an energy scale characterizing the strength 
of tunneling and A being an inverse length scale charac- 
terizing the degree of momentum non-conservation. We 
expect that the specific form of 7^ is unimportant as 
long as it includes the possibility of tunneling between 
different in-plane momentum states. Similar models that 
incorporate disorder in c-axis transport have been consid- 
ered by the authors of Refs. [23l25| who also find T 2 be- 
havior. In Ref. |25|. it is assumed that 7k depends only on 
the component of k parallel to the Fermi surface (imply- 
ing complete non-conservation of the perpendicular com- 
ponent). Under these assumptions, the leading temper- 
ature dependence of X~ (T) is quadratic, provided that 
the parallel component of k is conserved. It is however 
not easy to envision an interlayer scattering mechanism 
that would produce tunneling that is perfectly conserv- 
ing for the momentum component parallel to the Fermi 
surface while totally non-conserving in the perpendicular 
direction. A much more natural assumption, embodied 
by our choice Eq. (|15b|) . is to take 7k to be isotropic in 
the afr-plane. In the next section, we turn to the eval- 
uation of the c-axis penetration depth with this choice 
of tunneling matrix element. We will show that even 
isotropic tunneling produces a nearly T 2 dependence of 
S\ C (T) as a crossover behavior when the anisotropy of the 
Dirac spectrum near the node is taken into account. The 
same physics also produces the nearly x 2 doping depen- 
dence of A~ 2 (0), provided that we implement the charge 
renormalization factors as we did in the Sec. [H] 
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B. c-axis penetration depth 

The calculation of the c-axis penetration depth within 
the tunneling Hamiltonian formalism parallels that pre- 
sented in Appendix A for A a /,. The details can be found 
in Ref. [2f| and the result is given by 



1 

a! 



8e 2 ^T k 2 _ p T^F(k,c^(p,c), (16) 

k,p iuj 



As noted in Ref. [25|, this is equivalent to the standard 
result for the critical Josephson current through a weak 
link. This is sensible, since the way in which supercurrent 
is passed in the c-direction for weakly coupled layers is 
via the Josephson effect. 

We now augment Eq. I|16(l in the same manner as in 
Sec. [HI by making the replacement 



(17) 



k, P 



k, P 



With the choice made for Z k in Eq. © this again 
amounts to restricting the summations over k and p to 
within the vicinity of the nodes. Evaluating the sum over 
Matsubara frequencies, we obtain 



1 

A? 



2e 2 d^T k 2 _ p Z k Z 1 

k,p 



A k A p 
P E^Ep 



tanh \I3E^ + tanh \f3E p 
Ek + E p 
tanh |/3 i? k — tanh \f3E v 
E^-Ep 



(18) 



This equation contains a four dimensional momentum 
integral and one needs to employ numerical methods to 
obtain the full T and E c dependences of A c . We shall 
present such numerical results shortly. In order to elu- 
cidate the physics contained in this expression we first 
study the leading behavior analytically using scaling ar- 
guments in the limit of low T and low doping x [entering 
via E c oc x as in Sec. |nj)] when the physics is dominated 
by the nodal regions. 



C. c-axis penetration depth: scaling analysis 

Consider the low-T behavior of A 2 for T < E c < A . 
We shall denote the T-dependent correction to the stiff- 
ness by 5 A c (T) = A f T 2 (0)— A^ 2 (T). An explicit expression 
for SX C (T) may be obtained from Eq. I|18|) by replacing 
each tanh function by 1 — tanh (with the same argument). 
In the low-T limit, these functions restrict the momen- 
tum integrals to the vicinity of the nodes, so that we 
can take Z k ~ Zq everywhere. Linearizing near the four 
nodes and rescaling to the natural momentum variables 
as in Eq. (|8a|) we find that the tunneling matrix element 
Tj 2 becomes anisotropic, as illustrated in Fig. EI I n these 




FIG. 3: Schematic plot of the constant energy contours in 
momentum space in the vicinity of a node. On the left, con- 
tours at energy _E k satisfy _E k = Vpkf + v\k^ (i.e., they are 
ellipses). The tunneling matrix element Eq. I15b>t conserves 
momentum within a range A represented by the dashed cir- 
cle. By rescaling the plot so that the axes are vyk\ and v&ki, 
the constant energy contours are circles, but the circle rep- 
resenting the degree of momentum conservation has become 
distorted, indicating that the k\ component of the quasiparti- 
cle momentum is effectively conserved to a lesser degree than 
k 2 . 



natural "nodal" variables the momentum component per- 
pendicular to the Fermi surface (vpki) is conserved to a 
lesser degree than the parallel component (wa^)- This 
gives rise to three distinct scaling regimes depending on 
the magnitude of T with respect to the two natural scales 
of the problem, vpA and v A A. In what follows, we shall 
assume that v A <C Vp, as is known to be the case for 
cuprates>i& 

(i) For v A A <C vpA <C T thermally excited quasi- 
particles tunnel between the layers essentially as if mo- 
mentum were conserved, i.e., A is a small scale in this 
regime. In this limit we may approximate 1^ by a 
delta function and Eq. (|18f) becomes the same as the 
result Eq. J3J for A a fc, implying linear temperature de- 
pendence, S\ C (T) ~ T. 

(ii) For v A A <C T <C v f A the form of suggests 
that the dominant tunneling is characterized by v A k2 be- 
ing essentially conserved but v F k\ being essentially un- 
restricted. This is exactly the situation envisioned in 
Ref. [25|. leading to the T 2 behavior, but in our case it 
emerges as a crossover phenomenon. 

(iii) For T <C v A A <C vpA, extending this simple argu- 
ment (so that neither component is conserved) yields a 
vanishing result for SX C (T) due to the sign-change in the 
d-wave pair potential. Using a more careful Sommerfeld- 
like expansion detailed in Appendix [51 we find that 
5A C (T) oc T 3 in this regime. 

Summarizing, we have 



SXc(T) 



T 3 for T < v A A < v F A; 
T 2 for v A A «T« vpA; 
T for v A A < vpA < T. 



(19) 



The c-axis penetration depth data of Ref. [9j|10| exhibit 
power-law behavior consistent with T 3 crossing over to 
T 2 . The first two regimes of Eq. I|19|) indicate that, by a 
suitable choice of the parameter A, it may be possible to 
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fit this data with the present model. Before attempting 
such a fit, we turn to the E c dependence of X~ 2 (Q). 

Next, we carry out a similar analysis for A~ 2 (0), which 
has the explicit form of Eq. (|18fl but with each tanh func- 
tion replaced by unity. In analyzing SX C (T), we made use 
of the interplay between the way in which integrals were 
cut off by T and the way momenta were effectively con- 
served by ^_ p - For Aj (0), the momentum integrals are 
cut off by E c instead of T (via the factor Z^Z p ), but oth- 
erwise the previous analysis remains largely valid. In par- 
ticular, the intermediate (v A A < £ c < vpA) and high 
(v A A <C v-pA <C -E c ) E c regimes are analogous to their 
counterparts in Eq. i|19|) . At low E c , the naive analysis 
fails (as it did for 5X C (T) at low T); a more careful anal- 
ysis (described in Appendix lB|) shows that A~ 2 (0) cx E^ 
for E c — ► 0. We thus have 



El for E c < v A A < v F A; 
A~ 2 (0) ~ I El for v A A < E c < v F A; 

E c for v A A < v F A < E c . 



(20) 



Our results (|19|l and i|20|) indicate that the incoherent 
tunneling model with effective charge renormalization © 
qualitatively reproduces the trends in the c-axis penetra- 
tion depth summarized in Eq. @. The near-quadratic 
behavior in both temperature and doping observed in 
experimentSiifl emerges in our model as a crossover phe- 
nomenon involving the energy scales v A A <C vpA. We 
shall see below that the full numerical integration of Eq. 
(|18|l indeed reproduces the scaling behavior indicated in 
Eqs. I|19|) . I|20|) and furthermore gives excellent quantita- 
tive agreement with the experimental data. 

From our analysis above it should also be clear that 
if Zk fell to a small value proportional to y/x outside 
the nodal patch (instead of vanishing there) then X~ 2 (0) 
would pick up a contribution linear in x in all the regimes 
described in Eq. J5DJ, in disagreement with the data. 



D. c-axis penetration depth: numerical evaluation 

To facilitate numerical evaluation it is useful to rewrite 
Eq. (|18fl by converting the sums to integrations in the 
usual way and switching to relative and center-of-mass 
variables q = (k — p)/2 and Q = (k + p)/2. We have 



1 

a! 



16e 2 d 



d 2 Q 



{v F v A ) 2 J (2ttTi) 2 J {2irh) 



d 2 q 



"Q+q^Q-q 



(21) 



X T 2 

2q IQ 



l 



q||Q q|q Q 
'Q + q| 



|Q — q| tanh 



2T 



|Q + q| tanh 



IQ 



2T 



where we have also linearized in the vicinity of the nodes 
and rescaled the momenta v F Q\ — ► Qi and v A Q2 —> Qi 
and similarly for q. Such a rescaling effectively makes 
the argument of the matrix element factor anisotropic: 
In Eq. I|21|) its argument is q = (qi/v F , q2/v A ). 



1^,11) 
30 r -I 




r|=4.00 

— ■ ri=16.0 
. . . r]=50.0 



FIG. 4: Plot of h{T,n) for rj = 4, 16,50. To emphasize the 
crossover behavior in the power-law of I\ (r, n), the inset plots 
the associated logarithmic derivative a(j) = dln/i/dlnr. 



Before proceeding, we re-emphasize that incoherence 
in the tunneling matrix element is essential for obtain- 
ing the correct temperature dependence of A C (T) (and, 
with our choice for Z^, the correct doping dependence 
of A C (T)). For the case A — * 0, in which the spatially 
modulated tunneling matrix element becomes uniform 
on average, Eq. I|21|l yields a T-linear finite temperature 
correction X~ 2 (0) — X~ 2 (T). Thus, incoherence in the 
tunneling matrix element is crucial in what follows. 

It is convenient to express all quantities as dimension- 
ful prefactors multiplying dimensionless functions of di- 
mensionless parameters. We express (5A C (T) in this man- 
ner as 



-,??), (22a) 



oX c (l ) - — - Ji( , 

TT^V F V A hr y/VFV A A 

h(T,n)=4T 3 [ ePqe- 4r 'tih-& T i)fl(q), (22b) 



fi(q) 



1 



d 2 Q- 



ql-Ql 



1 



Q + q||Q q| q Q 
Q + qk 



|Q - q|(l - tanh 



|Q + q|(l - tanh 



IQ 



(22c) 



For simplicity in the above integrals we have set Z\^=Zq, 
an approximation valid as long as T <C E c . The integrals 
remain convergent due to the thermal Fermi factors. 

In the main part of Fig. 0] we display 1\ (r, rf) for 
T] = 4, 16, 50; other values of r\ display similar behav- 
ior. To focus on the power-law behavior of I(r,r)), in 
the inset of Fig. 21 we plot the logarithmic derivative 
q(t) = dhili/dhiT as a function of lnr, again for 
n = 4, 16, 50. The virtue of such a plot is that con- 
stant behavior of the logarithmic derivative at a particu- 
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V 2 (0) - K 2 w 2 ) 




FIG. 5: Plot of h{e Q ,ri) for r) = 4, 16,50. To emphasize the 
crossover behavior in the power-law of In(£ c ,r]), the inset plots 
the associated logarithmic derivative a(e c ) = dlnl2/d\ne c . 



lar value reflects power-law behavior of I(t, rj) with that 
value as the exponent. By examining this plot, there is 
clear evidence for r-linear behavior for r 3> ^/rj (in fact, 
the exponent approaches unity for much lower values of 
r) and r 3 behavior for r -C l/y 7 ??- At intermediate r 
the exponent is close to 2 although the expected plateau 
starts developing only for large anisotropy rj. We may 
thus conclude that <5A C (T) indeed exhibits behavior ex- 
pected from the scaling analysis presented above. 

We now analyze the doping dependence. It is conve- 
nient to express A,7 2 (0) in a form reminiscent of Eq. (|22a|l : 
i.e., as a product of a dimensionful prefactor multiplying 
a dimcnsionless function of suitably chosen dimensionless 
quantities. We thus have 



1 16e 2 Ad Zlt\ 

-h{- 



rj), (23a) 



(23b) 



A 2 (0) TTy/vFVA h 4 ^/wfwaA 
J 1 kp k+p 

vp- £ cPi-Pi) 2 Af-C=2-P2) 2 r;] 



where the subscript 1 on the integrations in Eq. 123 bl) 
indicate that the integration range is the unit disk, cor- 
responding to the discontinuous jump in our choice for 
Zk. In practice, for numerical convenience we shall re- 
place this hard cutoff with a Gaussian soft cutoff when 
performing numerical integrals. This corresponds to the 
choice Zy = Zq exp(— E^/ E'l) which leads to the same 
qualitative afe-plane behavior as discussed in Sec. [H] 

To ascertain whether the three regimes indicated in 
Eq. (|20|l are indeed realized, in the main part of Fig. [S] 
we display a numerical plot of ^(^c, v) f° r V = 4, 16, 50. 
The inset of Fig. we plot the logarithmic derivative of 
this quantity. For large e c , ^(^c, v) clearly exhibits linear 
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FIG. 6: Plot of fits to experimental data of Ref. |9J, 
ITc|| . The diamonds are \~ 2 (0) — A~ 2 (T) for various 
dopings having experimental T c values (top to bottom) 
20.2, 19.5, 18.2, 17.8, 16.4, 15. IK. The solid curve is our best 
fit using the parameters HA^ 1 = 120A and t± = 26meV. 
The inset is the same plot on a logarithmic scale, showing 
the changing power law of the experimental data and of our 
theoretical curve. 



behavior. For intermediate e c , there is clearly a regime 
of quadratic behavior that is more pronounced for the 
cases 77 = 16 and rj = 50 and the exponent is seen to 
approach 5 for low e c . The essential feature of hi^cV) 
is that it exhibits power-law behavior with an exponent 
near 2 for a considerable range of e c . The hope is that 
this quadratic dependence can, with a suitable choice of 
parameters, fit the known quadratic doping dependence 
of the c-axis penetration depth Eq. J5J. In the next sec- 
tion, we perform a detailed fit to c-axis dataAifi 



IV. c-AXIS PENETRATION DEPTH: DATA 
FITS 



In the present section, we attempt to fit the c-axis 
penetration depth data of Ref. |9lll0l| using the theoreti- 
cal formulas obtained in the preceding section. Thus, we 
assume that the experimental X~ 2 (T) ps A~ 2 (0) — SX C (T) 
with Aj 2 (0) and 5X C (T) being given by Eq. (|23a|l and 



Eq. 122a|) . respectively. Within this approximation to 
Eq. (|21|l . SX C (T) does not depend on E c , in agreement 
with the experimental observation that the finite-T cor- 
rection is universal (i.e. doping-independent) at low T. 
Thus, we begin by first fitting SX C (T) to the experimen- 
tal finite-T correction A,: 2 (0) - A^ 2 (T). 

The relevant fitting parameters are as follows: The pa- 
rameters t± and A characterize the way in which tunnel- 
ing occ urs in the c-axis direction; since the measurements 
of Ref. [H3 were done on one crystal, we shall take these 
to be global fitting parameters. The parameter Zq may 
be taken to be unity, as it only enters in Eq. (|23a|l and 
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T (K) 

FIG. 7: Plot of extracted values of the charge-renormalization 
parameter E c (diamonds) as a function of the experimental 
T c , showing a linear behavior as a function of doping level. 
The solid curve is a linear fit to these values, and has the 
form E c = 0.49T C /K + O.Ol(meV). 



Eq. in the product Z t±. We take v F = 1.8 eV Ail 

and d — 5.85A, although it is not known if vp changes 
for strongly underdoped samples. As discussed above, 
the way in which the power-law behavior of Eq. I|22a|) 
changes as quasiparticles are excited from the conden- 
sate depends on the magnitude of the anisotropy ratio 
7/ = tjf/^a; however, in practice we have not found the 
quality of the fits to be strongly dependent on this value. 
Thus, although the best fits (i.e. those minimizing the 
variance) were found with the value rj = 6.8, one cannot 
claim to extract a value for 77 using this technique. The 
converse of this is that although it is expected that 77 will 
vary somewhat for such underdoped samples (in a way 
that is not understood at present), the weak dependence 
of the quality of our fits on 77 validates our neglecting this 
doping dependence. 

In Fig. we show our best fits to the low-temperature 
values of A t T 2 (0) — \~ 2 (T) in the experiments of Rcf. 9, 
ITcj . The diamonds represent data curves for a partic- 
ular doping value. Each doping value is characterized 
by a particular T c which we take to be proportional to 
x. For clarity we have only displayed the highest doping 
values (i.e., T c = 20.2, 19.5, 18.2, 17.8, 16.4, 15. IK); the 
fits are equally good for lower doping values. The solid 
line is our best fit with the parameters h/A = 120 A and 
t± = 26meV. The fits work well at low T (despite the 
fact that the data is not a simple power law) but begins 
to deviate at high temperatures. (The roughly horizontal 
behavior occurs above T c ). We ascribe this discrepancy 
to fluctuation effects near T c in a given sample, as well 
as the fact that we have neglected the effect of on the 
finite temperature correction <5A C (T) in Eq. (|22al) . This 
restricts the validity of our calculations to low tempera- 
ture. As we saw in Sec. [H] for X ab , the deviations due to 
this effect are expected to become more pronounced at 
higher temperatures (see Fig- EJ - 



Having fit the temperature-dependent correction, the 
only remaining parameters are the values of E c corre- 
sponding to a particular doping. We extract these using 
Eq. (|23a|) for X~ 2 (0). In Sec. [H] we noted that to account 
for the a6-plane phenomenology, we must have E c cx x 
(and therefore E c cx T c ). To verify that this indeed holds, 
in Fig.[7|we plot (diamonds) the extracted best-fit values 
of E c for a given experimental T c from the Hosseini et 
al. results. The solid curve is a linear fit to these values, 
with the form E c = 0.49T C /K + O.Ol(meV). This agrees 
very well with the rough estimate of Sec. [H] where we 
found E c = 0.55T C meV/K from the ab plane data. This 
linear "Uemura"— relation is an important constraint on 
this theory and depicts the destruction of the Fermi sur- 
face as the Mott insulating phase is approached at low 
doping values. 

Finally, to illustrate the overall agreement of our model 
with the data, in Fig. [8] we plot the data of Hosseini et 
iilJLill for several representative doping values along with 
our curve fits. The agreement is strikingly good at low 
temperatures for all doping levels. We emphasize that 
all data sets are fit with a single set of parameters (listed 
in Fig. |SJ); the only parameter that varies is the cutoff 
energy according to E c — 0.49T c meV/K with T c being 
the actual measured critical temperature. 



V. CONCLUDING REMARKS 

In an effort to understand a6-plane and c-axis su- 
perfluid stiffness in underdoped cuprates we have con- 
structed a model incorporating incoherent tunneling be- 
tween the CuC>2 planes along with a strongly anisotropic 
charge renormalization factor Z^. Incoherent tunneling 
provides a mechanism for obtaining a non-linear temper- 
ature dependence of X~ 2 while the charge renormaliza- 
tion factor is introduced as means to model the funda- 
mental departure from the BCS theory in the underdoped 
regime where only doped holes contribute to the super- 
fluid. 

At this stage our choice for Z^ is purely phenomeno- 
logical and expands upon the suggestions of Refs. 
and [27|]. It is motivated by the observation^ that al- 
though the afr-plane and c-axis penetration depths are 
strongly doping dependent, the temperature-dependent 
corrections to these doping-dependent values are nearly 
doping- independent. At the coarsest level, this implies 
that quasiparticle effective charge Z^ ~ 1 for states near 
the nodes of the d-wave order parameter but is strongly 
reduced away from the nodes in a doping-dependent fash- 
ion. This doping dependence is governed by a cutoff en- 
ergy E c ~ x. By incorporating these two features into a 
model of the afr-plane and c-axis penetration depths, we 
were able to explain the low-temperature dat a 9 i 10 with a 
striking accuracy using a single set of input parameters 
and values of E c which vary remarkably linearly with 
doping level, as expected on the basis of the a6-plane 
phenomenology. 
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r(K) 

FIG. 8: Fit to data of Ref. [ij|[T(| (diamonds) using pa- 
rameters extracted in text. The T c values are T c = 
20.2, 18.2, 16.4, 12.1, 7.4K, (top to bottom) representing de- 
creasing effective doping. The parameters used are ?lA _1 = 
120A t ± = 26meV, n = v F /v A = 6.8 and E c = 
0.49T c meV/K. 

Taken together the above results lead to the notion of 
a "nodal protectorate" in which coherent BCS quasipar- 
ticles persist even as the system approaches the Mott- 
Hubbard insulating state near the half filling. This 
nodal protectorate is schematically depicted in Fig. 
which in addition illustrates how the protected regions 
shrink to essentially nothing as i — > 0. The exis- 
tence of this nodal protectorate is in addition supported 
by heat conduction 1 ^, scanning tunneling spectroscopy 36 
and photoemission experiments^. 

The existence of the nodal protectorate imposes a num- 
ber of stringent constraints on any microscopic theory 
describing the the underdoped regime. In particular any 
such theory must explain what protects the low-energy 
nodal excitations from the strong interactions that other- 
wise drive the electrons in the remainder of the Brillouin 
zone inert to applied electromagnetic fields. In addition 
one would like to understand what is the significance of 
the energy scale E c ~ x implied by our results. It is well 
known that RVB-type theories^ 9 - naturally explain the x- 
linear dependence of Pg h {0) but predict a coefficient of 
the T- linear term to go as x 2 , in strong disagreement 
with the experimental data. It has been claimed that 
the SU(2) version of the theory rectifies this problem 14 
but the physics of this is somewhat opaque and the sta- 
tus of this result remains unclear^ In addition it is not 
obvious that one would obtain the required c-axis behav- 
ior from this model. Approaches which seek to describe 
the Mott physics via Gutzwiller projection techniques ap- 
plied to the BCS wavefunction22i2i also obtain the correct 
Ps b {0) ~ x behavior. It is more difficult to address the 
finite-T properties within these models and the behav- 
ior of <5A a fc(T) is not known at present. However, these 
being real space techniques, it is not easy to envision a 
mechanism that would protect the nodal regions from 



the strong correlations imposed by the projection and 
one would naively expect that the result will suffer from 
the same problem as the RVB theories. 

It has been argued that the x and T dependence of p^ b 
can be explained within the Fermi liquid theory for su- 
perconductors with anisotropic Fermi surfac e) 39 ' 40 assum- 
ing a particular form of the interaction . We have 
checked that a most straightforward extension of this 
interaction does not reproduce the c-axis phenomenol- 
ogy; one would have to impose another set of constraints 
for the interlayer interactions to get the observed behav- 
ior which makes this route somewhat implausible in our 
view. One could imagine a competing order^i*^ in the 
particle-hole channel gapping out the Fermi surface away 
from the nodes, thus reducing the number of electrons 
participating in the superfluid in a manner consistent 
with the observed phenomenology^ii The advantage of 
this scenario is that nodal quasiparticlcs would automat- 
ically remain protected and the c-axis phenomenology 
would presumably also follow. On the other hand such 
competing orders necessarily break various space-time 
symmetries of the underlying system and such symmetry 
breaking should be easily observable if the competing or- 
der was strong enough to modify the superfluid response 
in accord with experiment. On balance we feel that the 
available evidence does not support this scenario^ 

In the phase fluctuation scenarios^^^ the low super- 
fluid density p s (0) ~ x enters as a phenomenological in- 
put. An appealing feature of these models is that in the 
superconducting state the phase fluctuations are gapped 
(with the gap of the order of p s ~ E c ) and thus do not 
affect the low energy nodal quasiparticles which remain 
sharp even as x — > 0. On the other hand the phase 
fluctuation models require some microscopic theory to 
describe the Mott physics that is ultimately responsible 
for the assumed low superfluid density. 

It would thus appear that none of the microscopic the- 
ories that exist at the present time satisfies all of the 
constraints implied by the phenomenological model ad- 
vocated in this paper. We may conclude that the physics 
of a <i-wave superconductor on the brink of becoming a 
Mott-Hubbard insulator remains an intellectual challenge 
awaiting future solution. 
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APPENDIX A: CALCULATION OF X ab 

In the present section, we compute the in-plane pene- 
tration depth for a BCS d-wave superconductor described 
by the Hamiltonian H in Eq. J3J with H mt — for now. 
This is related to the imaginary part of the in-plane con- 
ductivity (72 (oj) 



,35 



-i- = lim uja 2 (uj) 



(Al) 



The first step is to identify how the system couples to an 
applied electric field E x in the x direction. For the case 
of H having only nearest-neighbor hopping, the electro- 
magnetic current operator has the form 



tea 



r,cr v-\-dx , (7 



h.c), 



(A2) 



with a being the lattice spacing. The coupling to the 
electromagnetic vector potential A x (related to E x via 
A x = — ^E x e~ luJt ) enters by making the Peierls replace- 
ment t — > t exp(ieaA x /c) in Eq. ljA2fl . Expanding to 
to leading order in A, we have j x — j% + j x with the 
paramagnetic and diamagnetic currents being given, re- 
spectively, by 

Jx = ieat ^2( c l*c r +ax,<j -h.c), (A3) 

(J 

e 2 a 2 t 



fx = 



c 



C r,<j c r +ax. 



a +h.c.)A x . (A4) 



current is already linear in the electric field so that the as- 
sociated contribution to the conductivity can be directly 
read off by taking the expectation value of Eq. (|A6|I , 



D 



; E 



d 2 ek , t \ 

dkl ^k,^ C k,<r/' 



2^ EE 



%uj k 



dkl 



G(k,Lu) 



(A8) 
(A9) 



For the paramagnetic current, the usual leading-order 
perturbation theory result^ naturally leads to consider- 
ing the current-current correlator 



n(i/) 



dre- itJr (T T f(r)f(0)), 



2e2T EE(|r) [G(k,c)G(k,c-,) 

iuj k ^ x * 

+F(k,uj)F(k,iu~v)}, (All) 



(A10) 



where the in-plane Matsubara Green functions G(k, to) 
and FQc, ui) are given by the usual expressions 



G(k,w) 



F(k,w) = >iV , 



Ak 



2 ' 



(A12) 
(A13) 



and Ek = y / e 2 .+ A£. Performing the required Matsub- 
ara sums and combining the preceding expressions with 
Eq. ijAlf) . we find the penetration depth 



Our task is to find the conductivity <tl(w) for a single 
layer, which is defined by (j x (q,u>)) = cr(q,uj)E x (q,uj). 
The conductivity of the layered cuprate system will then 
be <j(uj) = ^ol(w) with n being the number of Cu02 
planes per layer and d being the interlayer spacing. The 
angle brackets represent the equilibrium average with re- 
spect to H, and henceforth we shall be concerned only 
with the q — * limit. In this limit, the Fourier trans- 
formed paramagnetic and diamagnetic current operators 
may be written as 



P — \ " ^ £k t 

Jx — e Qfc C k,cr C k,cr' 

k , (j x 

. d _ ?e 2 \ - d 2 e k t 
Jx ~ u ^ dkl '' k T '' k 

k,o- * 



(A5) 
(A6) 



where we have displayed the expression for general ek 
(and will continue to do so henceforth). The frequency 
dependent conductivity per layer is 



<j l (0,lo) = -[D + TP(uj)], 



(A7) 



1 e 2 n v - 



d 2 e k 
dkl 



1 - ^- tanh l(3E k 



2 

en 
~2d~ 



^E 



9ek 

dk x 



2 1 

sech -pE\c, 



(A14) 



with (3 = 1/T the inverse temperature. This last expres- 
sion agrees with the results for X a b in the lattice model of 
a d-wave superconductor obtained by other authors i^^ii 
As written, the diamagnetic term in Eq. I|A14|I has 
contributions from the entire Brillouin zone while the 
paramagnetic term is dominated by the nodal regions 
(Ek — > 0). For our purposes it will be convenient to re- 
cast the former into a form where it is also dominated 
by the nodal regions. This can be accomplished by in- 
tegrating by parts in the first term of Eq. i|A14(l . Wc 
obtain 



1 

A ab 



2 

en 



k L 



dfk 

dk x 

1 

Ek" 



<9e k <9A k A k e k 



dE k 



tanh 2^-^k- 



(A15) 



where D and II represent the diamagnetic and param- 
agnetic parts of the response. As usual, the diamagnetic 



This expression is mathematically equivalent to 
Eq. (|A14|) and gives X a b as a fc-space sum dominated 
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by the nodal regions. In addition, Eq. I|A15(1 has the 
desirable property of explicitly yielding zero superfluid 
stiffness in the normal limit, At = 0. Finally we notice 
that the second term in the first line of Eq. (|A15|) 
is generally small for a <i-wave gap. In particular it 
vanishes identically in the nodal approximation that we 
employ in Sec. [HJ Thus, we shall ignore this term. 



APPENDIX B: ASYMPTOTIC BEHAVIOR OF 

8X C {T) AND A~ 2 (0) 

In the present section we sketch the derivation of the 
low T behavior of 5X C (T) and the low E c behavior of 
A~ 2 (0), i.e., the first lines of Eq. JTHJ and Eq. (JUJ, 
respectively. For we use the smooth form = 
Z exp(-El/E%) introduced in Sec. lIIIDI We will omit 
unimportant prefactors and, for simplicity, set v-p — v& — 
1. It will be clear from the derivation that relaxing this 
last assumption will not change the results. We begin 
with SX C (T), which is given by 



SX C (T) cx 



d 2 p d 2 k _ (k _ p) 2 /A 2 k 2 p 2 
{2irh) 2 (27rS)2 e kp 
p(l - tanh(fc/2T)) 

p2 _ f.2 



(Bl) 



Henceforth, the polar coordinate expressions of k and 
p are (k,9) and (p,4>), respectively. Our strategy is to 
utilize the following Sommerfeld expansion: 



dkf{k)(l -tanhfc/2T) 



(B2) 



dk[f(0) + fc/'(0) + •••](!- tanhfc/2T), 



which relies on the sharpness of the function 1 — 
tanh/c/2T near k = for T — ► and is valid provided 
that the function f(k) is sufficiently smooth near k = 0. 
In the present case, the function / is given by (again, up 
to numerical prefactors) 



/(*)= I dpdH0 kp2 t 6 ^ ^ P)2/A2 

7} — K 



(B3) 



The first thing to note is that, due to the factor of k 
in Eq. I|B3(1 (in turn arising from the k in the measure 
of Eq. CE)), /(0) = 0. In addition, /'(0) = 0. The 
easiest way to see this is to note that, among the terms 



in f'(k), the only one that can possibly be nonzero at 
k = is the one for which the derivative acts on the k in 
the numerator of the fraction in Eq. l|A13jl . Thus, 



no) 



dp d(f> d6 sin 9 sin <f> e 



-P 2 /A 2 



0. 



(B4) 



as the angular integrals vanish by symmetry. Clearly, we 
must consider /"(0). A direct evaluation of /"(0) does 
indeed yield a finite result, so that the leading contribu- 
tion comes from the third term in Eq. (|B2|) . Evaluating 
the integral over k in this term then yields SX C (T) oc T 3 
for T — > yielding the first line of Eq. lj*H?|l . 

Next, we turn to the e c dependence of A~ 2 (0). It is 
simplest to consider the scaled function l2(s c ,r)), which 
may be written as 



/ a (Ec,l) = I d 2 k I d?p k -^^- k2 /^-V 2 /^ 



kp k + p 



xe 



-(k- P ) 2 



(B5) 



Equation l|B5|) is obtained from Eq. (|23b(l by replacing 
the integrals over the unit disk with a smooth Gaussian 
cutoff and then rescaling k — > k/e c (and similarly for p). 
We have also set rj = 1 for convenience. We proceed in 
the same manner as for the case of SX C (T) above, writing 
I 2 (e c , 1) as 



9(k) 



dkg(k)e 



(B6) 



dk[g(0) + kg'(0)+---}e- k2 ^, (B7) 



p dp d(f> d9 



k sin 9 sin i 

k + p 
,p-(k- P ) 2 



-p 2 /e 



(B8) 



Owing to the fact that the Gaussian function is sharply 
peaked near k = 0, we have in Eq. (|B7p approximately 
evaluated Eq. (|B6p by Taylor expanding g{k) in anal- 
ogy with Eq. I|B2|) . As in the previous case, g(0) = 
and <7'(0) = 0, requiring the evaluation of g"(0). Thus, 
the calculation proceeds exactly as above, with one im- 
portant difference: The function g(k) contains a factor 
q-p / e c j n the integrand. It is straightforward to verify 
that the integration over p then gives g"(0) cx e\. Com- 
bining this with a factor of e 3 arising from the third term 
in Eq. I|B7(I . we have that ^(^o 1) oc for e c — ► 0, im- 
mediately giving the first line of Eq. (|2Tj|l . 
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